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Coherent large scale jets that are not forced directly at the jet scale are a prominent 
feature of rotating turbulence. These jets arise and are supported by systematic orga- 
nization of the turbulent Reynolds stresses. Understanding the mechanism producing 
the required eddy momentum flux convergence, and how the jets and associated eddy 
field mutually adjust to maintain a steady jet structure at finite amplitude, constitute 
fundamental theoretical problems. Stochastic Structural Stability Theory (SSST) pro- 
vides a framework for understanding the emergence and equilibration of jets based on 
a statistical mean state dynamics of turbulence that is closed at second order. SSST 
reveals a new manifold of turbulent equilibria and linearization of the SSST dynamics 
about these equilibria allows formulation of a structural stability theory that predicts the 
emergence of jets from homogeneous turbulence. In this work we test predictions of this 
statistical mean state dynamics for both formation and equilibration of turbulent jets us- 
ing quasi-linear and nonlinear simulations of barotropic turbulence on a beta-plane. We 
find that SSST accurately predicts both the initial emergence and the finite amplitude 
equilibration of jets in beta-plane turbulence. 
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1. Introduction 

Spatially and temporally coherent jets are a common feature of turbulent flows in 
planetary atmospheres with the banded winds of the giant planets constituting a familiar 
example (Vasavada & Showman 2005). Fjortoft (1953) noted that the essentially two 
dimensional large scale turbulence in planetary atmospheres, which results from the 
constraints of rapid rotation and strong stratification, would tend, given approximate 
conservation of energy and vorticity, to transfer turbulent energy to the largest scale 
available. This argument provides a conceptual basis for understanding the observed 
tendency for formation of large scale structure from small scale turbulence in planetary 
atmospheres. However, the observed large scale structure primarily takes the specific 
form of zonal jets and, moreover, the scale of these jets is distinct from the largest scale 
in the flow. Rhines (1975) argued that the observed jet scale results from the arrest of the 
inverse cascade at the length scale, y/u//3, where j3 = d//dy is the meridional gradient 
of planetary vorticity and u is the root mean square velocity in the turbulent fluid. 
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In Rhines' interpretation this scale is identified with the scale at which the turbulent 
cascade is intercepted by the formation of propagating Rossby waves. More recently it 
was proposed that jets in beta-plane turbulence arise from the tendency of turbulence 
to mix potential vorticity so that it becomes homogenized in limited regions forming 
potential vorticity staircases. The risers of these staircases correspond to thin prograde 
jets located at the latitudes of steep potential vorticity gradients separating parabolic 
retrograde jets corresponding to the well mixed steps of the staircase (Baldwin et al. 2007; 
Dritchel & Mclntyre 2008). Another explanation for jet emergence has been advanced 
by Miller (1990) and Robert & Sommeria (1991). They associate formation of large 
scale jets with the tendency of two dimensional turbulent flows to form structures of 
maximum entropy while respecting energy and vorticity conservation. These statistical 
mechanical considerations predict the most likely structures to be zonal jets (Bouchet 
& Venaillc 2012). Emergence of large scale structure in planetary flows has also been 
related to the existence of fixed points of the nonlinear planetary scale flow dynamics. 
Fixed point solutions were initially advanced to explain large scale persistent phenomena 
in the atmosphere such as blocking (Charney & DeVore 1979; Legras & Ghil 1985). These 
fixed point structures are solutions to the large scale laminar flow dynamics and do not 
include the influence on the dynamics of the interaction between the large scale structure 
and the smaller scale turbulence. Yet another proposed explanation for jet formation is 
modulational instability in which a zonal flow arises from nonlinear interaction of discrete 
Rossby waves (Gill 1974; Connaughton et al. 2010). 

SSST provides an alternative explanation for jet formation that is fundamentally based 
on the interaction between jets and their associated field of turbulent eddies (Farrell & 
Ioannou 2003). In this theory, jets initially arise from a linear instability of this interaction 
and finite amplitude jets are identified as the nonlinear extension of these instabilities. 
SSST combines the full dynamics of the zonal mean flow with the second order statistics 
of the turbulent field obtained from a stochastic turbulence model (STM). The STM has 
been shown to provide accurate eddy statistics for the midlatitude atmosphere (Farrell 
& Ioannou 1993, 1994, 1995; Zhang & Held 1999; O'Gorman & Schneider 2007). These 
coupled equations determine the evolution of the mean zonal flow under the influence of 
its consistent field of eddies. The SSST equations predict jet formation from instabilities 
of the stationary solutions of this coupled dynamics and the instability of the stationary 
solutions as a function of parameters determines the bifurcation structure of the jet for- 
mation process. Srinivasan & Young (2012) developed a continuous formulation of SSST 
and used it to obtain analytic expressions for the emerging SSST instability in homo- 
geneous beta-plane turbulence. A related quasi-linear approach to understanding finite 
amplitude turbulent equilibria in planetary atmospheres is the second order cumulant 
truncation method (CE2) (Marston et al. 2008; Marston 2012). 

The quasi-linear (QL) approximation to the full nonlinear dynamics (NL) results when 
the perturbation-perturbation interactions are parameterized in the perturbation equa- 
tions while interaction between the perturbations and the zonal mean flow is retained 
fully in the zonal mean equation. The SSST is based on QL with the additional assump- 
tion that the perturbation dynamics is that of the STM, which is equivalent to an infinite 
ensemble of perturbations evolving under QL. Because the SSST equations provide an 
autonomous and fluctuation-free dynamics of the statistical mean turbulent state, these 
equations facilitate analytic study of turbulent equilibria and their stability. 

Previous studies demonstrated that unstable jets maintained by mean flow forcing 
can be equilibrated accurately using QL dynamics (Schoeberl & Lindzen 1984; DclSolc 
& Farrell 1996; Marston et al. 2008). In contrast to these studies, in this work we in- 
vestigate the emergence and equilibration of jets from homogeneous turbulence without 
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coherent external forcing. SSST predicts that infinitesimal perturbations with zonal jet 
form organize homogeneous turbulence to produce systematic upgradient fluxes giving 
rise to exponential jet growth and eventually to the establishment of finite amplitude 
equilibrium jets. Specifically, the SSST equations predict initial formation of jets by the 
most unstable eigenmode of the linearized SSST dynamics. The formation of zonal jets 
from initially homogeneous beta-plane turbulence was studied with the methods of SSST 
by Farrell & Ioannou (2007); Bakas & Ioannou (2011, 2012) and by Srinivasan & Young 
(2012). In agreement with SSST, Srinivasan & Young (2012) found that NL simulations 
exhibit jet emergence from a homogeneous turbulent state with subsequent establishment 
of finite amplitude jets. In this paper we use the dynamics of the full NL equations and 
their QL approximation to examine further the mechanism of emergence and equilibra- 
tion of jets from turbulence, concentrating on the effects of perturbation-perturbation 
nonlinearity on jet bifurcation and equilibration and on the influence of perturbations in 
exciting the manifold of SSST modes with jet structure. 



2. Formulation of nonlinear barotropic dynamics on a beta-plane 

Consider a beta-plane with x and y Cartesian coordinates along the zonal and the 
meridional direction respectively. The nondivergent zonal and meridional velocity fields 
are expressed in terms of a streamf unction, tp, as u — —d y ip and v = d x ip. The vorticity 
on the beta-plane is q + 2£l + f3y, where q = d x v — d y u is the relative vorticity and 2S} + /3y 
is the planetary vorticity. NL is governed by the beta-plane barotropic vorticity equation: 

d t q + ud x q + vd y q + 0v = y/fF -rq-P 4 A 2 q . (2.1) 

The flow is dissipated with linear damping at rate r and hyperviscosity with coefficient v^. 
Periodic boundary conditions are imposed on x and y with periodicity 2nL. In the above 
equation, distances have been nondimensionalized by a chosen length scale L and time by 
T = L/U, where U is a chosen velocity scale. We choose L — 5000 km and U — 40 ms _1 
for which the nondimensional time unit is T = 1.5 day and the nondimensional terrestrial 
midlatude value of j3 is 10. Turbulence is maintained by stochastic forcing. The structure 
of the imposed forcing is specified by the function F and its amplitude by the scalar /. 
This spatially and temporally homogeneous forcing will be specified later. 

Velocity fields are decomposed into zonal-mean components (indicated upper case) and 
perturbations (indicated lowercase and primes) so that the zonal velocity is U(y,t) + 
u'(x,y,t), the meridional velocity is v'(x,y,t), and the vorticity is Q(y,t) + q' (x,y,t). 
Averaging in the zonal direction is denoted with a bar. Taking the zonal average of 
Eq. (2.1) we obtain two coupled equations for the evolution of the mean flow, U, and the 
associated perturbation field: 

d t U = W-r m U , (2.2a) 

d t q' = -Ud x q' + {U yy - P) d x ^ - rq 1 - ^ A 2 a' + F e + y/fF , (2.26) 

where 

F e = (9, (Wq 7 ) - d y (v'q'/j - d x (u'q') (2.3) 

is the nonlinear term representing the perturbation-perturbation interactions. In the 
above equation ip' = A" 1 ^', where A -1 is the inverse of the Laplacian A = d xx + dy y . 
The mean flow and perturbation Eqs. (2.2) define the nonlinear system (which will be 
referred to as NL). In Eqs. (2.2) the damping rates r and r m may differ. Hyperviscosity 
is included in the perturbation equations for numerical stability and forcing is applied 
only to the perturbation equations. 
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3. Formulation of quasi-linear barotropic dynamics on a beta-plane 

The QL approximation of NL (2.2) is obtained by parameterizing the perturbation- 
perturbation nonlinearity in Eq. (2.26) as a state independent stochastic excitation and 
retaining the nonlinear Reynolds stress forcing, v'q', in the mean zonal flow equation to 
obtain the quasi-linear (QL) system: 

d t U = W-r m U , (3.1a) 

d t q' = -U8 x q' + (U yy - /3) d x i>' - rq' - U 4 A 2 q' + y/fF , (3.16) 

where F gives the spatial and temporal structure of the stochastic excitation. 

The perturbation-perturbation nonlinearity, F e , redistributes energy and enstrophy 
among the eddy scales. It can be shown that both the complete system (2.2) and its QL 
approximation (3.1) conserve both energy and enstrophy in the absence of forcing and 
dissipation. 

Because in Eqs. (3.1) U is not a function of x, and the perturbation-perturbation 
nonlinearity has been neglected, the perturbation vorticity, q' , can be represented as a 
sum of non-interacting zonal harmonics 



q.'{x,y,t) = Re 



E 

U-l 



(3.2) 



where Re denotes the real part and 1, . . . , Nk are the zonal wavenumbers excited by the 
stochastic forcing in (3.16). The structure of the stochastic forcing is expanded as: 



F(x,y,t) = Re 



E 

lk=l 



ikx 



(3.3) 



The stochastic forcing Fk(y,t) — J2t=i ^kt(y) €hi(t) excites each of the Fke(y) (to be 
specified) with coefficient £ke(t)- The £^ are temporally (5-correlated and independently 
excite both the zonal harmonics and the meridional structures and are prescribed by: 



&<(*)) = , (^(*)Cn(*0) = SkrnSin 5(t - t') 



where the angle bracket denotes an ensemble average over realizations. 
The two-point perturbation vorticity covariance is defined as: 

C(x a ,y a ,x b ,y b ,t) = (q'(x a: y a ,t)q'(x b ,y b: t) 



(3.4) 



(3.5) 



Because of the independence of the excitations, and the separability of the linear 
dynamics (3.16) in x, the covariance is homogeneous in the zonal direction at all times 
and can be expanded in the difference coordinate, x a — x b , as: 



C(x a - y a , y b , t) = ^ Re 



Re 



E 

,fc=i 

E 

,fc=i 



qk(y a ,t)qt(y b ,t) 



ik(x a —Xb) 



C k (y a ,y b ,t) e 1 *^-**) 



(3.6) 



where Ck{y a , y b , t) is the latitudinal covariance of the vorticity field at zonal wavenumber 
k and time t. 

The spatial covariance of the forcing is defined as 



Q(x a ,x c ,y b ,y d ) = ^F(x a ,y b ,t)F(x c ,y d ,t)^ , 



(3.7) 
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Q(x a - x c , y b - y d ) = Re 



N k N v 

Y,Y.^ e ik(x a - Xa )+i(y b -y d )} 
k=l 1=1 



1=1 



(3.8) 



At zonal wavenumber k, the corresponding latitudinal spatial covariance is: 

N v 

Qk(ya,yb) = J^Qkie u ^-yK (3.9) 



4. Formulation of SSST barotropic dynamics on a beta-plane 

Taking the Fourier transform of Eq. (3.16) and discretizing the fields on a meridional 
grid we obtain the corresponding matrix form: 

d t q k = A k (U)q k + ^fF k , (4.1) 

with 

A k (U) = -ik [U - {U yy - /3I) A^ 1 ] -rl- v 4 A 2 . (4.2) 

The state, q k ; the zonal mean flow, U; and the Fourier coefficient of the forcing, F k ; 
are column vectors with elements the values at the N y discretization points, yj. In the 
matrix A k , A k = D 2 — k 2 l, with D 2 the discretized d yy operator, / the identity matrix, 
A k l the inverse of A k , U the diagonal matrix with diagonal elements U and U yy the 
diagonal matrix with diagonal elements D 2 U. With this discretization the latitudinal 

covariance Ck(y a ,yb,t) of Eq. (3.6) becomes the matrix C k = (q k q k f, with f denoting 
the Hermitian transpose. 

The meridional spatial structure of the forcing, F k , is specified by a matrix F kl so that 
for each k, the white noise ^ k i excites the spatial structure given by the £-th column of 
Ffc. Consequently, the forcing of the k-t\i wavenumber is written as: 



■ k\ j 



Ny 

E 

1=1 



\F k ] f M (t) 

L J jl 



(4.3) 



The structure of the forcing is determined by F kl and the forcing structures have been 
normalized, as will be described. 

The dynamical equation for the ensemble average of the covariance of the perturbations 
in Eq. (4.1), C k , satisfies the time dependent Lyapunov equation: 

d t C k - A k (U)C k + C k A k (U)i + fQ k , (4.4) 

as in Farrell & Ioannou (1996). The overall forcing amplitude is controlled by the pa- 
rameter / and the spatial covariance of the forcing enters this equation as Q k = F k F k , 
or equivalently in terms of the spectral coefficients of the forcing covariance, Q k £, as: 



bii 



N y 

£< 

e=i 



hi e 



'^{Vb-Vd) 



(4.5) 



The covariances, C k , determine the second order statistics of all perturbation fields. 
For example, the ensemble average of the vorticity flux, 

' %&)) , (4-6) 



(v'q 1 ) 



k=l 
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is expressed in terms of the covariances, C k , as: 



vq=~J2o vecd Im (\Ak l 9k9h 



k=l 

N k 



o vecd 



k=l 



(4.7) 



where is the column vector of the values of (v'q'} and Im denotes the imaginary part. 
We denote with vecd (M) the column vector consisting of the diagonal elements of matrix 
M. 

Under the ergodic assumption that the zonal average of the perturbation covariance 
can be obtained from the ensemble average, we substitute v'q' = (v'q') in Eq. (3.1a) to 
obtain an equation for the evolution of the zonal mean, which is linear in the C k : 

N k 



d t [/ = -^|vecd lm(A^C k 



fc=i 



r m U . 



(4.8) 



The system of equations formed by (4.4) and (4.8), define the SSST dynamics for the evo- 
lution of the zonally averaged flow, U, and its associated ensemble average perturbation 
covariances, C k : 



k=l 



r m U , 



d t C k = A k (U)C k + C k A k {U)t + fQ k 



(4.9a) 
(4.96) 



5. Specification of the stochastic forcing 

Three stochastic forcing structures will be used in our investigation of the correspon- 
dence of SSST, QL and NL dynamics. 

The first independently excites a set of zonal wavenumbers. This forcing was first 
used by Williams (1978) to parameterize excitation of barotropic dynamics by baroclinic 
instabilities. This method of forcing was also used by DelSolc (2001) in his study of upper- 
level tropospheric jet dynamics and in the study of jet formation using SSST dynamics 
by Farrell & Ioannou (2003, 2007) and Bakas & Ioannou (2011). This stochastic forcing 
is spatially homogeneous but not isotropic. 

The second is an isotropic ring forcing concentrated near a single total wavenumber. 
This forcing structure has been used extensively in studies of beta-plane turbulence (cf. 
Vallis & Maltrud 1993) and was also been used in the recent study of Srinivasan & Young 
(2012). It was introduced by Lilly (1969), in order to isolate the inverse cascade from the 
forcing in a study of two dimensional turbulence. 

The third is an isotropic ring forcing, similar to the second, but with the forcing 
distributed over a wider annular region around the central wavenumber. 

Recall that the stochastic forcing, 



F(xi,yj,t) = Re 



fc=i i=\ 



£ki(t)e 



ikx , 



(5.1) 



is correlated in y by the columns of the matrix F k . For the non-isotropic forcing (NIF) 
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= Cfe 


' e -(y,-ye) 2 /(2s 2 ) + e -(y 3 -2K-y e ) 2 /{2s 2 ) 


e -( yj +2Tr-ye) 2 /(2s 2 ) 




it 







this meridional structure of the forcing correlation, in a channel periodic and of length 
2n in y, is specified by: 



(5.2) 



Each forcing is Gaussian correlated in y, with half-width s = 0.2/\/2 and we force zonal 
waves k = 1, . . . , 14. Because the stochastic forcing is <5-correlated in time, the energy 
injection rate does not depend on the state and can be independently specified (cf. 
Appendix B). The normalization constant, Cfe, in (5.2) is chosen so that each k is excited 
equally and one unit of energy is injected in total, so that / is the total energy injection 
rate in the NL, QL or SSST simulations. 
The isotropic ring forcing is specified by: 



F k =cw kl {K)e u yi 



(5.3) 



with K = y/W 
integers and 
specified by 



t 2 . Because the channel is square with length 2ir, both k and I are 
integers and only integer wavenumbers are excited. The first isotropic forcing (IRFh) is 



\K 
\K 



K, 



^Sk f 
>5k f 



(5.4) 



For IRFh we choose Kf — 14 and 6k f = 1. The second forcing (IRFg) is specified as 



w k i = c exp 



2 5kj 



(5.5) 



For IRFg we choose Kf = 14 and 5k f — 8. The normalization constant, c, is chosen so 
that the total energy injection rate is unity (cf. Appendix B). IRFh and IRFg are both 
spatially homogeneous and nearly isotropic in a finite doubly periodic domain. They 
approach exact isotropy as the domain size increases. 

Plots of the corresponding power spectra together with instantaneous realizations of 
the stochastic forcing both in vorticity and streamfunction are shown in figure 1 for the 
three types of forcing structures. The thin IRFh ring forcing primarily excites vortices of 
scale 1/Kf in the flow that are evident in both the vorticity and streamfunction fields, 
while IRFg produces a streamfunction field dominated by large scale structure. 



6. Stability of the SSST homogeneous equilibrium state 

Because of the homogeneity of the stochastic forcing, vq = in the absence of symme- 
try breaking by a mean flow, and the SSST equations admit the homogeneous equilibrium 
solution: 

N k , 

U E = , C £ = ^Cf with Qf = |-Q fe , (6.1) 

k=l 

for — 0. The stability of this equilibrium state, (U E ', C s ), is determined by the linear 
dynamics of the evolution of perturbations (SU,6Ci, . . . ,5C^ k ) to this equilibrium. These 



8 



N. C. Constantinou, P. J. Ioannou and B. F. Farrell 




Figure 1. Stochastic forcing spatial wavenumber spectra and examples of realizations in physical 
variables. Shown for NIF with zonal wavenumbers k = 1, . . . , 14 and s = 0.2/V2 are: (a) 
enstrophy spectrum, Q k e, (b) vorticity, (c) streamfunction. For IRFh at Kf = 14 and Skf — 1: 
(d) enstrophy spectrum Qke, (e) vorticity, (f) streamfunction. For IRFg at Kf = 14 and 6k f = 8: 
(g) enstrophy spectrum Qke, (h) vorticity, (i) streamfunction. 



perturbations satisfy the linear equations: 



N k 



d t SU = - ~ v ecd ImiA^SQ 



k=l 



r m SU 



d t SC k = Af5C k + 5C k {Af ) f + SAkCi + C% (SA k ) f , 



(6.2a) 
(6.26) 



which are obtained by linearizing the SSST system (4.9) about (U E , C E ). In the above 
equations A E = A k (U E ) is the operator linearized about the equilibrium mean zonal 
flow, U E , and SA k is SA k = — ifc \5U — (SU) yy A^ 1 ] ■ Eigenfunctions of the form, 



(6U,5C 1 ,...,6C Nh )e° 



(6.3) 
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satisfy the following eigenvalue problem: 



a5U = -^]-vecd Im (A^SCk) 



k=l 



r m SU 



a5C k = A E 6C k + 5C k (A E ) f + 5A k C E + Cjf (SA k ) f 



(6.4a) 
(6.46) 



The equilibrium (U E , C E ) is stable when all the eigenvalues of system (6.4) satisfy ay — 
Re (a) ^ 0. Otherwise, the SSST equilibrium is unstable and the SSST system bifurcates 
to a structurally distinct attractor. 

The eigenfunctions and eigenvalues of Eqs. (6.4) are readily obtained in the case of the 
homogeneous equilibrium state (6.1). Because of the homogeneity of this equilibrium, 
the mean flow part of the eigenfunctions are sinusoidal with form SU n — sin (ny). This 
perturbation of the mean flow consists of n pairs of alternating prograde and retrograde 
jets. For the case (3 = the eigenvalues can be obtained analytically (Bakas & Ioannou 
2011), otherwise both the eigenfunction, SC k ^ n , and the eigenvalue, o~(n), can be obtained 
iteratively by introducing SU n = sin (ny), which specifies SA ki . n in (6.46), obtaining 5C k , n 
from the resulting Sylvester equation: 



A E - a(n)l] SC k . n + SC k . n (A E ) f = - \5A k>n C E + C E (5A k J 



(6.5) 



introducing 6C k , n in Eq. (6.4a) and continuing this iteration until convergence is obtained. 

We now investigate the stability of the homogeneous equilibrium state (6.1) as a func- 
tion of the forcing amplitude, /. In the absence of forcing, / = 0, Eqs. (6.4) decouple and 
all mean flow perturbations decay with cr(n) = —r m . In the presence of forcing the mean 
zonal flow perturbations are coupled to the perturbation covariance and this coupling 
may give rise to an SSST instability. We will demonstrate that such a structural insta- 
bility occurs for forcing amplitudes / > f c resulting to transition to an inhomogeneous 
equilibrium with finite amplitude zonal jets. 

Growth rates, o>, are shown in figure 2 for NIF as a function of the meridional 
wavenumber of the mean zonal flow perturbation, n, and of the forcing amplitude, /, and 
similarly for IRFh in figure 3. For both types of forcing, instability occurs for f > f c over 
a band of mean flow wavenumbers, n. In all cases the o~(n) associated with the largest 
growth rate has zero imaginary part. 

Growth rate curves are shown in figure 4 for various forcings and damping amplitudes. 
In the next section the predictions of this SSST stability analysis will be compared with 
the corresponding QL and NL simulations. While QL and NL simulations can not provide 
theoretical predictions for the bifurcation structure of turbulent equilibria, we wish to 
examine the circumstances under which the underlying bifurcation structure revealed by 
SSST is reflected in the QL and NL realizations. 



7. Bifurcations predicted by SSST and their reflection in QL and NL 
simulations 

We examine the reflection of the SSST structural instability in NL and QL simulations 
by comparing the evolution of the domain averaged energy of the zonal flow: 

E m (t) = J ^-JJ 1 -U 2 dxAy. (7.1) 

The amplitude of the zonal flow is measured with the zonal mean flow index (zmf ) defined 
as zmf = E m /(E m + E p ), where E m is the time averaged energy of the zonal mean flow, 
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NIF k = 1, 



14 r = 0.1 
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Figure 2. Growth rate, o>, of the SSST eigenfunction with jet structure SU n = sin (ny) as a 
function of meridional wavenumber, n, and forcing amplitude, /, for NIF with k — 1, . . . , 14. The 
stability boundary (oy = 0) is marked with thick solid line and the unstable region is shaded. 
The least forcing required for instability, f c — 8.16 x 10 -5 , occurs at n = 6. Other parameters 
are: r = 0.1, r m = 0.01, = 10. 




Figure 3. Growth rate, oy, of the SSST eigenfunction with jet structure SU n = sin (ny) as 
a function of meridional wavenumber, n, and forcing amplitude, /, for IRFh at Kf — 14 with 
5k f — 1. The stability boundary (oy = 0) is marked with thick solid line and the unstable region 
is shaded. The least required for instability, f c = 5.82 x 10 , occurs at n — 8. Other parameters 
are: r = 0.1, r m = 0.01, P = 10. 



given in Eq. (7.1), and E p is the time average of the domain average kinetic energy of 
the perturbations 

E p (t) ^^-JJ 1 - ( u ' 2 + z/ 2 ) dxdy . (7.2) 
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/ = 1.5/„, r = 0.1 / = 20/ c , r = 0.1 




/ = 2.5f c , r = 0.1 / = 10/ o r = 0.01 




Figure 4. Growth rate, er r , of the SSST eigenfunction with jet structure 8U n = sin (ny) as a 
function of jet meridional wavenumber, n. For NIF with k — 1, . . . , 14 shown are the growth 
rates for (a) / = 1.5/ c and r = 0.1, (b) for / = 20/ c and r = 0.1. For IRFh at K f = 14 with 
Skf — 1 shown are: (c) the growth rates for / = 2.5/ c and r = 0.1, (d) for / = 10/ c and r = 0.01. 
Other parameters are: r m = 0.01, = 10. 



7.1. Correspondece of jet emergence in NL, QL and SSST for moderate damping 

The case of NIF with r = 0.1, r m = 0.01 and / = 1.5/ c is plotted in figure 5. For this 
forcing and parameters, SSST predicts that maximum growth occurs for the sinusoidal 
zonal mean flow perturbation with n — 6, as seen in figure 4a. Introduction of this 
eigenfunction in the SSST dynamics produces initially the predicted exponential growth, 
followed by approach to an equilibrium consisting of 6 finite amplitude jets. This equi- 
libration demonstrates that the SSST dynamics includes the mechanism of stabilization 
of the SSST instability. Corresponding simulations with the QL and NL dynamics using 
the same parameters for the damping and the forcing reveal nearly identical jet growth 
followed by finite amplitude equilibration. A Hovmoller diagram showing the develop- 
ment of the mean flow associated with the SSST instability is shown in figure 6, together 
with corresponding diagrams for the development of the associated mean flow in the QL 
and NL simulations. Although no theoretical prediction of this bifurcation behavior can 
be made directly from the QL or the NL dynamics, they both reveal the bifurcation 
structure obtained from the SSST analysis. 

At larger forcing amplitudes the SSST jet formation instability typically equilibrates at 
a wavenumber smaller than that of the instability. An example for NIF at supercriticality 
/ = 20/c is shown in figure 7. While the jets emerge initially with zonal wavenumber 
n = 6 (cf. figure 4b), in agreement with the prediction of SSST instability, the zonal 
flow eventually equilibrates at wavenumber n — 4, following a series of mergers that are 
evident in the Hovmoller diagrams. This behavior can be rationalized by noting that 
as the jet amplitude increases at a fixed wavenumber, violation of the Rayleigh-Kuo 
stability criterion occurs, at which point, by transitioning to a lower wavenumber, the 
flow forestalls occurrence of an inflectional instability (Farrell & Ioannou 2007). However, 
SSST stability analysis of the finite amplitude equilibria near the point of jet merger 
reveals that these mergers coincide with the inception of a structural instability, which 
precedes the occurrence of hydrodynamic instability of the jet (Farrell & Ioannou 2003). 

A plot of the zonal mean flow indices (zmf ) in the NL and QL simulations and in the 




Figure 5. Development of mean flow energy, E m , for NL and QL simulations starting from a 
state of rest with NIF at k = 1, . . . , 14 and / = l-5/ c . Also shown is the exponential growth and 
subsequent equilibration of E m under SSST dynamics with initial condition the most unstable 
SSST eigenfunction. This figure demonstrates that SSST dynamics accurately predicts nonlinear 
jet equilibration at finite amplitude as well as the initial instability and exponential growth of 
the jets. Parameters are /3 = 10, r = 0.1, r m = 0.01. 



associated SSST equilibria are shown in figure 8 as a function of forcing amplitude for 
NIF forcing of k = 1, . . . , 14 with damping parameters r = 0.1, r m = 0.01. Consistent 
with the specific examples shown in figure 6 and figure 7, this figure demonstrates that 
the SSST bifurcation structure is paralleled by the QL and NL simulations over a wide 
range of forcing amplitudes. It also demonstrates that both the NL and QL simulations 
equilibrate at finite amplitude with nearly the same structure as the SSST equilibria. 

Similar agreement is obtained with IRFh at Kf — 14 with Skf = 1 and damping 
parameters r = 0.1, r m = 0.01 for which the case with / = 2.5/ c is shown in figure 9. 
Maximum growth occurs at zonal wavenumber n = 8 in this case (cf. figure 4c) and 
again the SSST predictions for both initiation of the instability and its saturation are 
reflected in the QL and NL simulations. The bifurcation diagram for this IRFh forcing and 
damping is plotted in figure 10. For this case also the critical amplitude, f c , is correctly 
predicted by SSST stability analysis, but for supercritical forcing amplitudes the zonal 
flow predicted by SSST and QL is stronger and the perturbation energy smaller than the 
corresponding zonal flow and perturbation energy obtained in the NL simulations. This 
difference between the NL and both the SSST predictions and QL simulations at high 
levels of forcing, corresponding to high superctiticality, is in part accounted for by the 
spreading of the spectrum in NL as shown in figure 10, where the zmf index that results 
when the SSST is forced by the NL spectrum at / = 20/ c is shown. However, the major 
part of this difference is due to the nonlinear perturbation-perturbation interactions in 
NL that disrupt the upgradicnt momentum transfer. This disruption is accentuated by 
the peculiar efficiency with which the narrow ring forcing, IRFh, gives rise to vortices as 
can seen in figure 1 (middle panel) . We verify that this narrow ring forcing is responsible 
for depressing NL equilibrium jet strength at high superctiticality by broadening the 
forcing distribution in the form of IRFg for which Kf = 14 and Skf = 8. Using this 
modified forcing we obtain agreement between SSST, QL and NL simulations with the 
same parameters as those used in figure 10, as is shown in figure 11. 
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Figure 6. Hovmoller diagrams of jet emergence in the NL, QL and SSST simulations shown in 
figure 5. Shown for the NL, QL and SSST simulations are U(y,t) (left panels) and for the NL 
and QL simulations the vorticity fields at t = 2000 (top and middle right panels). Also shown 
are the equilibrium jets in the NL (dash-dot-black), QL (dashed - red), and SSST (solid - blue) 
simulation (bottom right panel). This figure shows that SSST predicts the structure, growth 
and equilibration of weakly forced jets in both the QL and NL simulations (cf. figure 4a). 



7.2. Correspondence of jet emergence in NL, QL and SSST for weak damping 

For sufficiently small perturbation damping, r, the SSST critical forcing amplitude, f c , for 
jet emergence differs from the critical forcing amplitude, /i NL ' ) , at which finite amplitude 
and persistent zonal mean flows form in the NL simulations. For example, for perturbation 
damping r = 0.01 the NL simulations bifurcate at /J NL ^ ~ H/ c under NIF forcing with 
k = 1, . . . , 14 and at /i NL ^ ~ 4/ c under IRFh forcing at Kf = 14 and 5kf = 1. Similar 
behavior was noted by Srinivasan & Young (2012). The bifurcation diagram for the NIF 
forcing is shown in figure 12 and for the IRFh forcing in figure 13. The reason for the 
difference in the NL and SSST bifurcation curves will be explained in section 8. Typical 
structure of the emerging jets for the IRFh forcing is shown in figure 14. 
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Figure 7. Hovmoller diagrams of jet emergence in NL, QL and SSST simulations with NIF 
forcing of k — 1, . . . , 14 with / = 20/ c . Shown for the NL, QL and SSST simulations are U (y, t) 
(left panels) and for the NL and QL simulations the vorticity fields at t = 2000 (top and middle 
right panels). Also shown are the equilibrium jets in the NL (dash-dot-black), QL (dashed - 
red), and SSST (solid - blue) simulation (bottom right panel). This figure shows that SSST 
predicts the structure, growth and equilibration of immoderately forced jets in both the QL and 
NL simulations (cf. figure 4b). Parameters are r — 0.1, r m = 0.01 and /3 = 10. 



8. Development of the turbulence spectra in nonlinear simulations 
and their influence on SSST stability 

Both QL and SSST dynamics exclude interactions among perturbations and include 
only the non-local interactions between zonally averaged quantities, with k — 0, and per- 
turbations, with k 0. Therefore, there is no enstrophy or energy cascade in wavenum- 
ber space in QL and SSST dynamics and the equilibrium state in SSST, (U E — and 
C E = f Sfc=i 0fe/2r), has a spectrum proportional to the spectrum of the forcing, which 
is fully prescribed by the forcing covariances, Qk- However, this is not true in NL dynam- 
ics which includes perturbation-perturbation interactions producing enstrophy/energy 
cascades. For example, in NL an isotropic ring forcing is distorted as time progresses 
by the nonlinear perturbation-perturbation interactions becoming concentrated at lower 
wavenumbers and forming the characteristic dumbbell shape seen in beta-plane turbu- 
lence simulations (Vallis & Maltrud 1993), as shown in figure 15. It is to be expected 
that the statistical equilibrium state obtained by NL dynamics will be different from the 
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Figure 8. Bifurcation structure comparison for jet formation in SSST, QL, and NL. Shown is 
the zmf index of jet equilibria for NIF at k = 1, . . . , 14 as a function of the forcing amplitude 
f/fc for the NL simulation (dash-dot and circles), the QL simulation (dashed and dots) and the 
SSST (solid and diamonds). Parameters are /3 = 10, r = 0.1, r m = 0.01. 

corresponding SSST equilibrium state due to this difference in the spectrum. The equi- 
librium spectrum obtained by taking a time average of long NL simulations, in which the 
k = component of the flow has been set to zero, is shown in figure 16 for various forcing 
amplitudes and damping rates. The departure of the NL spectra from the spectra of the 
QL and SSST equilibria is evident and this difference depends on the amplitude of the 
forcing and the linear damping coefficient, r. A necessary condition for obtaining agree- 
ment between NL and both SSST and QL dynamics is that the equilibrium spectrum of 
the SSST and QL dynamics be close to the equilibrium spectrum in the NL dynamics so 
that the stability analysis is performed on similar states. 

The equilibrium spectrum obtained in the NL simulation does not differ appreciably 
from the spectrum of the forcing for IRFh at Kj = 14 with 5k f = 1 and linear damping 
r = 0.1 (cf. figure 16). In such cases SSST successfully predicts the critical parameters 
for the emergence of jet structure in NL simulations as shown by the bifurcation diagram 
in figure 10 and, for the case of NIF, in figure 8. 

SSST dynamics can be shown to accurately predict the emergence of zonal mean flows 
in the more weakly dissipated regime with r — 0.01, for which the NL spectra depart 
significantly from the spectra of the SSST homogeneous equilibrium, by performing the 
SSST stability analysis on the equilibrium spectrum obtained in the NL dynamics prior 
to emergence of the zonal mean jets. For example, with IRFh at Kf — 14 and 5k f = 1 
and perturbation dissipation r = 0.01, formation of persistent finite amplitude zonal 
mean flows occurs in the NL simulations in the range 3/ c < f < 4/ c instead of at f c 
(cf. figure 13). We perform an SSST stability analysis with forcing covariance associated 
with the spectrum fQki = 2r{\qkt\ ), where (|<to| 2 ) is the enstrophy spectrum obtained 
from the NL simulations in which the k = component was removed. We find that the 
equilibrium state with this spectrum is more stable than the original isotropic ring SSST 
equilibrium. This can be seen in figure 17 in which the growth rate of the SSST instability 
is shown as a function of the meridional wavenumber of the zonal mean perturbations 
for various forcing amplitudes. This calculation reveals that the NL-equivalent SSST 
equilibrium is marginally stable for / = 2.8/ c , and the NL modified spectrum produces 
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Figure 9. Hovmoller diagrams of jet emergence in NL, QL and SSST simulations with IRFh 
forcing at / = 2.5/ c and with Kf = 14, Skf = 1. Shown for the NL, QL and SSST simulations 
are U(y,t) (left panels) and for the NL and QL simulations the vorticity fields at t — 1500 (top 
and middle right panels). Also shown are the equilibrium jets in the NL (dash-dot- black), QL 
(dashed - red), and SSST (solid - blue) simulation (bottom right panel). This figure shows that 
SSST predicts the structure, growth and equilibration of jets obtained in both the QL and NL 
simulations (cf. figure 4c). Parameters are r = 0.1, r m = 0.01 and ft = 10. 



SSST instability for / ^ 2.8/ c . This result agrees with the NL simulations which show the 
emergence of persistent zonal mean flows for / ^ 2.8/ c . Moreover, SSST stability analysis 
on this NL-equivalent equilibrium state predicts that the zonal mean flow instability will 
occur at meridional wavenumbcr n = 5 or n = 6, while the maximal instability for the 
isotropic ring spectrum that characterizes the corresponding SSST equilibrium occurs 
at the higher wavenumber n — 8 (cf. figure 4d). In agreement with these predictions, 
zonal mean flows emerge in the NL simulations at / = 10/ c with n = 6 while in both 
the QL and SSST simulation zonal mean flows emerge with n = 9, and equilibrate 
to a finite amplitude zonal mean flow with n = 6, as shown in figure 14. We find in 
addition that the reduced amplitude of the zonal flows seen in the NL simulations (cf. 
figures 12 and 13) can also be traced to the modification of the forcing spectrum by the 
nonlinear processes which reduces the forcing from the eddy field. The zmf that results 
in SSST using the NL modified spectrum is shown in figures 12 and 13 to agree with 
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Figure 10. Bifurcation structure comparison for jet formation in SSST, QL, and NL. Shown 
is the zmf index of jet equilibria for IRFh with Kf = 14 and 5k f = 1 as a function of the 
forcing amplitude f/f c for the NL simulation (dash-dot and circles), the QL simulation (dashed 
and dots) and the SSST (solid and diamonds) for IRFh (SSSTa), as well as the SSST for 
the equivalent spectrum to that obtained in NL as a consequence of modification of the IRFh 
spectrum by the nonlinear cascade (SSSTb) with forcing amplitude / = 20/ c . This figure shows 
that the bifurcation to jet formation in NL and QL simulations is captured by SSST and that 
SSST overestimates the NL jet amplitude when the forcing is high, even when account is taken 
for the nonlinear modification of the spectrum. Parameters are /3 = 10, r — 0.1, r m = 0.01. 




Figure 11. Bifurcation structure comparison for jet formation in SSST, QL, and NL. Shown is 
the zmf index of jet equilibria for IRFg with Kf = 14 and 5k f — 8 as a function of the forcing 
amplitude f/f c for the NL simulation (dash-dot and circles), the QL simulation (dashed and 
dots) and the SSST (solid and diamonds). Parameters are /3 = 10, r = 0.1, r m = 0.01. 



zmf obtained in the NL simulations f. This demonstrates that differences between the 

f The NL modified spectrum used to obtain the SSSTb bifurcation curve shown in figure 12 
was obtained from the NL equilibrium spectrum at / = 4/ c . The spectrum used in figure 13 was 
obtained by similar means using an NL simulation at / = 2/ c (both with the k=0 component 
removed). The results are not sensitive to the choice of spectrum as long the broadening of the 
spectrum has been represented. 
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Figure 12. Bifurcation structure comparison for jet formation in SSST, QL, and NL. Shown is 
the zmf index of jet equilibria for NIF at k = 1, . . . , 14 as a function of the forcing amplitude 
f/fc for the NL simulation (dash-dot and circles), the QL simulation (dashed and dots) and 
the SSST NIF (SSSTa), as well as the SSST for the equivalent spectrum to that obtained in NL 
as a consequence of modification of the NIF spectrum by the nonlinear cascade (SSSTb). The 
bifurcation in the NL simulations occurs at / C (NL) « ll/ e . Thisfi gure shows that the structure of 
the jet in NL simulations is captured by the SSST if account is taken of the nonlinear modification 
of the spectrum. Parameters are f3 — 10, r = 0.01, r m = 0.01. 
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Figure 13. Bifurcation structure comparison for jet formation in SSST, QL, and NL. Shown is 
the zmf index of jet equilibria for IRFh at Kf = 14 with 8k j — 1 as a function of the forcing 
amplitude f/f c for the NL simulation (dash-dot and circles), the QL simulation (dashed and 
dots) and the SSST IRFh (SSSTa), as well as the SSST for the equivalent spectrum to that 
obtained in NL as a consequence of modification of the IRFh spectrum by the nonlinear cascade 
(SSSTb). The bifurcation in the NL simulations occurs at / C (NL) « 4/ c . This fi gure shows that 
the structure of the jet in NL simulations is captured by the SSST if account is taken of the 
nonlinear modification of the spectrum. Parameters are /3 = 10, r — 0.01, r m = 0.01. 
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Figure 14. Hovmoller diagrams of jet emergence in NL, QL and SSST simulations with IRFh 
forcing at / = 10/ c and with Kj — 14, Skf = 1. Shown for the NL, QL and SSST simulations 
are U(y,t) (left panels) and for the NL and QL simulations the vorticity fields at t — 1500 (top 
and middle right panels). Also shown are the equilibrium jets in the NL (dash-dot- black), QL 
(dashed - red), and SSST (solid - blue) simulation (bottom right panel). This figure shows that 
SSST predicts the structure, growth and equilibration of jets obtained in both the QL and NL 
simulations (cf. figure 4d). The jet maintained in the NL simulation has an amplitude which is 
10% less than the amplitude of the jet of the QL simulation. Parameters are /? = 10, r = 0.01, 
r m = 0.01. 



predictions of SSST stability theory and NL simulations for the emergence of zonal mean 
flows under weak damping results from the distortion of the equilibrium spectrum by the 
nonlinear interactions in the NL simulations. 

We can appreciate the control that the perturbation spectrum exerts on jet emergence 
by examining the case of forcing amplitude / = 2/ c with IRFh. For this forcing ampli- 
tude, the SSST is structurally unstable and both SSST and QL simulations equilibrate 
to finite amplitude jets, as shown in figure 18. SSST stability analysis performed on the 
NL modified spectra (cf. figure 19) shows that for this forcing amplitude, the asymp- 
totic spectrum is subcritical to structural instability (cf. the SSST growth rates on the 
approximately asymptotic t = 250 spectrum in figure 20). Consistent with this asymp- 
totic stability, the NL simulation does not support coherent persistent jets (cf. figure 18). 
However, while in the NL simulation no coherent jets are seen, clear indications of jet 
structure can be seen, although these jets are intermittent in both space and time. In 
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Figure 15. Evolution toward equilibrium of the energy spectrum, Eke = \ \<jkt\ 2 /(k 2 + £ 2 ), 
in NL under IRFh forcing with Kf — 14 and 8kf = 1. At t = the spectrum is that of the 
homogeneous SSST equilibrium. Top panels: The spectrum at t = 1, 20, 37 for / = 0.75/ c . 
Bottom panels: The spectrum at t = 1, 20, 37 for / = 2/ c . The upscale energy cascade 
produces the familiar dumbbell spectrum. The distortion of the spectrum increases with /. The 
parameters are r — 0.01, r m = 0.01 and /3 = 10. 
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Figure 16. Equilibrium enstrophy spectrum, log (2r|g«| 2 ), of NL simulations, in which pertur- 
bation-perturbation interactions are included but the k = component is excluded for various 
damping rates. The example is for IRFh forcing at / = 2/ c and with Kf = 14, 8k / = 1. Shown 
is the spectrum: (a) for linear damping coefficient r — 1 (left panel), (b) r = 0.1 (middle panel) 
and r = 0.01 (right panel). The critical, f c , is a function of r and is obtained from SSST. All 
spectra have been normalized by their maximum value. This figure shows that for relatively 
small damping the perturbation spectra in NL simulations are close to the SSST perturbation 
spectra while for weak damping the equilibrium spectrum in NL differs substantially from that 
in SSST. Other parameters are /3 = 10 and r m = 0.01. 



order to understand these observations we perform SSST stability analysis on the spec- 
trum as it develops in the NL simulation. The developing spectrum, shown in figure 19, 
is obtained from a simulation initialized with a perturbation and mean flow free state 
and the corresponding SSST growth rates associated with the developing spectrum are 
shown in figure 20. By t = 20 the NL spectrum, having assumed the isotropic ring struc- 
ture of the forcing, becomes SSST unstable. This structural instability is evident in the 
incipient formation of jet structure in the NL simulation at t — 50 in figure 18. As the 
spectrum further evolves, the SSST growth rates fall and the spectrum supports only 
stable modes for t > 120, with decay rates asymptotically approaching those correspond- 
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Figure 17. SSST jet instability growth rates for the equilibrium spectrum obtained in the NL 
dynamics subjected to IRFh forcing at Kf = 14, Skf = 1 for r = 0.01 (shown in the right panel 
of figure 16). Shown are the growth rates, o>, as a function of the meridional wavenumber, n, 
with forcing amplitude / = 2/ c , / = 2.8/ c , / = 4/ c and / = 10/ c . Also shown are the growth 
rates for the SSST equilibrium with IRFh forcing at / = 2/ c (solid-star line). This figure shows 
that the nonlinear interactions have modified the original ring spectrum making it more stable 
so that the modified flow state becomes unstable at / = 2.8/ c and n — 6. With this modification, 
the NL simulations agree with the predictions of the SSST analysis as shown in figure 13. Other 
parameters are f3 — 10 and r m = 0.01. 



ing to t = 250 in figure 20. For these asymptotic perturbation spectra, SSST stability 
analysis reveals a manifold of decaying modes among which the least damped modes are 
jets with meridional wavenumber n = 5, 6, 7, 8. There is clear evidence of zonal structure 
in the NL simulations in the form of a superposition of the n = 5, 6, 7, 8 least damped 
modes of the SSST stability analysis (cf. figure 17). 

These damped modes are excited by the fluctuations associated with realizations of 
the stochastic forcing in NL simulations. Similar excitation of intermittent jets associ- 
ated with these damped modes arise in QL simulations with SSST subcritical stochastic 
forcing. An example simulation (cf. figure 21) shows these intermittent jets both in QL 
and NL simulations under SSST subcritical forcing. Confirmation that these intermittent 
jets result from fluctuation excitation of the SSST damped modes is given in the bottom 
panel of figure 21, where the intermittent jets resulting from stochastic forcing of the 
SSST modes themselves are shown. This diagram was obtained by plotting 



U(y,t)=Re 



N 

,71=1 



(8.1) 



with a n independent red noise processes, associated with the damping rates, —a(n), of 
the first N = 15 SSST least damped modes. These a n are obtained from the Langevin 
equation, 

^=o-{n)a n + t{t) , (8.2) 
at 

with £(t) a 5-correlated complex valued random variable. 

The fluctuation-free SSST simulations reveal persistent jet structure by the inception 
of the SSST instability only at supercritical forcing. However, in the QL and NL simula- 
tions fluctuations excite the damped manifold of modes predicted by the SSST analysis 



22 



N. C. Constantinou, P. J. Ioannou and B. F. Farrell 



NL U(y,t) , f/f c =2 



■ 0.02 
0.01 




q(x, y, t = 2000) 



I 



500 1000 1500 

QL U(y,t) , f/f c =2 



2000 



-0.01 
-0.02 




0.05 



g(x, y, i = 2000) 



0.05 




2 4 6 







0.05 c 






4 






a 






2 






-0.05 



-0.05 




0.05 



Figure 18. Hovmoller diagrams of jet emergence in NL, QL and SSST simulations with IRFh 
forcing at / = 2/ c and with Kf = 14, Skf = 1. Shown for the NL, QL and SSST simulations 
are U(y,t) (left panels) and for the NL and QL simulations the vorticity fields at t — 2000 (top 
and middle right panels). Also shown are the equilibrium jets in the NL (dash-dot- black), QL 
(dashed - red), and SSST (solid - blue) simulation (bottom right panel). At / = 2/ c the NL 
simulation is subcritical but the n — 6 structure of the least stable SSST mode is discernible 
(cf. figure 16). This figure shows that stable SSST eigenmodes can be excited by fluctuations in 
NL simulations in the subcritical regime. Parameters are /3 = 10, r = 0.01, r m = 0.01. 



to exist at subcritical forcing amplitudes. This observation confirms the reality of the 
wave/mean flow interaction manifold of stable modes. Moreover, the emergence of struc- 
ture in subcritical simulations may be related to the observation of the so-called latent 
zonal jets in the oceans (Berloff et al. 2009, 2011). 

In NL and QL simulations these stable modes predicted by SSST are increasingly 
excited as the critical bifurcation point in parameter space is approached, because their 
damping rate vanishes at the bifurcation. The associated increase in mean flow energy 
on approach to the bifurcation in NL and QL simulations, obscures the exact location of 
the bifurcation point in NL and QL simulations compared to the fluctuation-free SSST 
simulations for which the bifurcation is exactly coincident with the inception of instability 
(i.e. figures 8, 10 and 11). 

SSST dynamics not only predicts the emergence of zonal jets but also predicts finite 
amplitude the equilibrium jets that result after the initial jet formation instability. These 
finite amplitude jets correspond to fixed points of the SSST dynamics. Here we compare 
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Figure 19. Evolution of the ensemble average enstrophy spectrum, |<Jfcfj 2 , for a set of NL 
simulations with isotropic ring forcing at wavenumber Kj — 14 , 5kf — 1, and forcing amplitude 
/ = 2/ c . Parameters are /3 = 10, r = 0.01, r m = 0.01. 
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Figure 20. Instantaneous SSST growth rates, oy, as a function of jet meridional wavenumber, 
n, for equilibria with the enstrophy spectra shown in figure 19. The spectrum of the evolving 
NL simulation becomes SSST unstable at t w 20 and stabilizes again at t w 120. Parameters 
are /3 = 10, r = 0.01, r m = 0.01. 



the SSST predictions of finite amplitude jet structure with the results of NL and QL 
simulations. For the case of NIF (cf. figure 7), as well as for for IRFh and IRFg, the 
jet equilibria of the SSST dynamics are similar to the jets obtained in QL and NL 
simulations. An example for IRFh strongly forced with / = 100/ c and with damping 
r = 0.01 is shown in figure 22. This case demonstrates that a 20% difference in the zmf 
index alone between the NL, QL and SSST (cf. figure 13) is not reflected in a substantial 
disparity between jets in NL, QL and SSST simulations. 
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Figure 21. Hovmoller diagrams of intermittent jet structure in NL and QL simulations at sub- 
critical forcing / = 0.8/ c . Shown are U(y,t) resulting for: NL (top) and QL (middle) simulations 
and the U(y,t) that results from random excitation of the SSST damped modes (bottom panel). 
These plots were obtained using IRFh forcing with Kf = 14, Skf = 1 and r — 0.1 (cf. figure 10). 
This figure shows that the manifold of SSST damped modes are revealed by being excited in 
the fluctuating NL and QL simulations. Other parameters are j3 = 10, r m = 0.01. 




Figure 22. The mean jet structure, U{y), for NL (dash-dot), QL (dashed) and SSST (solid) 
simulations for IRFh at Kf = 14 with Skf — 1 for / = 100/ c and perturbation damping r = 0.01. 
There is good agreement between the jet structure in the NL simulation and the jets obtained 
with QL and SSST, despite the difference in the zmf index (cf. figure 13). The other parameters 
are /3 = 10, r m = 0.01. 
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9. Conclusions 

In this work we compared the predictions of SSST for jet formation in barotropic beta- 
plane turbulence with QL and NL simulations. We found that the bifurcation structure 
predicted by SSST for the emergence of zonal jets from a homogeneous turbulent state 
is confirmed by both QL and NL simulations. Moreover, we showed that the finite am- 
plitude equilibrium jets found in NL and QL simulations are as predicted by the fixed 
point solutions of SSST. Obtaining this agreement between NL and both SSST and 
QL simulations required that the modification of the turbulent spectrum caused by the 
perturbation-perturbation nonlinearity in NL be accounted for in the specification of the 
stochastic forcing in QL and SSST. These results confirm that jet emergence in barotropic 
beta-plane turbulence can be traced to the cooperative mean flow/perturbation instabil- 
ity that is captured by SSST. Moreover, the physical reality of the manifold of modes 
arising from cooperative interaction between incoherent turbulence and coherent jets, 
which is revealed by SSST, was verified in this work by relating the excitation of the 
spectrum of damped SSST modes to observations of latent jets in both QL and NL sim- 
ulations. SSST provides an autonomous, deterministic nonlinear closure of the turbulent 
dynamics at second order that provides a vehicle for further analytic investigation of the 
dynamics of coherent structures in turbulent flows. 
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Appendix A. Numerical details and parameters 

Periodicity in the square channel [0, 2n] x [0, 2tt] implies that the permissible zonal and 
meridional wavenumbers, (k,£), are integers and the NIF and IRFh or IRFg forcing is 
limited to excite only these integer wavenumbers. The coefficient of linear damping for 
the mean flow in all simulations is r m = 0.01. The perturbations are either damped at 
the rate r = 0.1 or at the rate r = 0.01. In all simulations f3 = 10. 

Both the nonlinear simulation (NL) Eq. (2.1) and the quasi-linear simulations (QL) of 
Eqs. (3.1) are done with a pseudospectral Fourier code. The maximum resolved wavenum- 
bers are fc max = N x /2 and £ max = N y /2. For the time integration we use a fourth order 
Runge-Kutta method (RK4) using a Godunov step for integrating the stochastic forcing. 

The coefficient of is chosen for numerical stability as 



where KjI = 



(N x /2f + (N y /2) 



vl = KiJ*' (A1) 

2 

is the maximum resolved total wavenumber in the 

fourth power and St is the time step. In all our calculations St = 2.5 x 10 -3 , which implies 
v A = 1.86 x 10~ 7 . 

The stochastic forcing at the j-th time step is obtained from: 

i km (t-) = '^fi , (A 2) 
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where rjj and 9j, j — 1,2,... are independent normally distributed random variables with 
zero mean and unit standard deviation. The phase of £ km is uniformly distributed on 
[0, 2tt]. The v8t in the denominator ensures that the forcing maintains constant variance 
as temporal resolution changes and the y2 ensures that (; km has unit variance, i.e. both 
ensure that in the limit St —> we obtain: J (£,km(t)^ km (t')) dt' = 1. 



Appendix B. Calculation of the energy injection rate 

The domain averaged energy injection rate by the stochastic forcing, F, is: 

dx dy 



L X Ly 



(Bl) 



We Fourier expand all fields as in (3.2). When y is discretized, the Fourier coefficient of 
the streamfunction, ipk(y, t), of the k zonal wavenumber, is related to vorticity coefficient, 
Qk(y,t), by V'fe(i) = A k l q k (t) where i/)fc and q k are column vectors and A^ 1 the inverse 
of the matrix representation of d yy — k 2 . The solution of Eq. (4.1) in matrix form is: 



q k {t) = exp (A k t) q k {0) + 7/ / exp \A k (t - s)l F k £ k {s) ds , 



(B2) 



where £ k is the column vector with elements [£ k ] e = £,ki- Assuming ergodicity, the zonal 
average in Eq. (B 1) can be replaced by an ensemble average and then converted in matrix 
form, obtaining: 

dy 



-77 



-4>'F 




fc=i 



NyA 



(B3) 



with Tr denoting the trace of a matrix. Because the noise processes are independent and 
temporally delta correlated, satisfying (£fcm(£)£/c'm'(0) = $kk'&mm'5(t — t'), we obtain, 
using (B2), that both terms of the sum in (B3) are independent of the state of the 

system and equal to — 1-\/7 Tr ^ / 7 1 F fc F^ . Consequently, energy is injected at the rate: 



k=l k=l 

in the NL, QL and SSST simulations. 



J_ 

47V,. 



Tr (^fe 1( ?fc) > 



(B4) 
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